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SUMMARY 

Coupled problems with various combinations of multiple physics, scales, and domains can be found 
in numerous areas of science and engineering. A key challenge in the formulation and implementation 
of corresponding coupled models is to facilitate communication of information across physics, scale, 
and domain interfaces, as well as between iterations of solvers used for response computations. In a 
probabilistic context, any information that is to be communicated between subproblems or iterations 
should be characterized by an appropriate probabilistic representation. In this work, we consider 
stochastic coupled problems whose subproblems involve only uncertainties that are statistically 
independent of one another; for these problems, we present a characterization of the exchanged 
information by using a reduced chaos expansion with random coefficients. This expansion provides 
a reduced-dimensional representation of the exchanged information, while maintaining segregation 
between sources of uncertainty that stem from different subproblems. Further, we present a measure 
transformation that allows stochastic expansion methods to exploit this dimension reduction to obtain 
an efficient solution of subproblems in a reduced-dimensional space. We show that owing to the 
uncertainty source segregation, requisite orthonormal polynomials and quadrature rules can be readily 
obtained by tensorization. Finally, the proposed methodology is demonstrated by applying it to a 
multiphysics problem in nuclear engineering. Copyright © 2000 John Wiley & Sons, Ltd. 
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1. Introduction 

In the fields of climate research, combustion, and renewable energy and in many other 
critical areas of science and engineering, models and simulations depend on a common 
base of mathematical formulations, algorithms, and implementations involving multiphysics, 
multiscale, and multidomain characteristics. Demand for predictive computational results in 
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these areas has motivated the development of uncertainty quantification (UQ) approaches for 
coupled problems with various combinations of multiple physics, scales, and domains. 

Probability theory provides a rigorous mathematical framework for UQ. The first step in a 
probabilistic UQ analysis typically involves the use of mathematical statistic;^ methods [1, 2] 
for the probabilistic characterization of the uncertain features associated with a model as 
one or more random variables, random fields, random matrices, or random operators. The 
second step is to use the model to map this probabilistic representation of inputs into a 
probabilistic representation of responses. This can be achieved in several ways, for example, 
by using Monte Carlo sampling techniques [3] and stochastic expansion methods. The latter 
typically involve the representation of the predictions as a polynomial chaos (PC) expansion. 
Several approaches are available to calculate the coefficients in this expansion, for example, 
embedded projection [4, 5], nonintrusive projection [5], and collocation [6-10]. 

A key challenge in the formulation and implementation of coupled models is to facilitate 

communication of information across physics, scale, and domain interfaces, as well as between 
iterations of solvers used for response computations. This information can comprise physical 
properties, energies, solution patches, and other quantities. Although the number of sources of 
uncertainty can be expected to be large in most stochastic coupled problems, we believe that 
exchanged probabilistic information often resides in a considerably lower dimensional space 
than the sources themselves. In stochastic multiphysics and stochastic multiscale problems, 
the exchanged information can be expected to be of low effective stochastic dimension when 
it consists of a solution field that is smoothed by a forward operator and when it is obtained 
by summarizing fine-scale quantities into coarse-scale representations, respectively. 

In a previous paper [11], we had thus proposed the use of a dimension-reduction technique, 
namely, the Karhunen-Locve (KL) decomposition, to represent exchanged probabilistic 
information. In a subsequent paper [12], we had presented a measure-transformation technique 
that allows stochastic expansion methods to exploit this dimension reduction to obtain a 
computationally efficient solution of subproblems in a reduced-dimensional space. In [11] 
and [12], we had considered a general setting that allowed the uncertainties in the data of 
the various subproblems to be statistically dependent on one another. 

In this paper, we present an extension of our previous work. This extension is applicable 
to stochastic coupled problems whose subproblems incorporate only sources of uncertainty 
that are statistically independent of one another; for these stochastic coupled problems, we 
present a characterization of the exchanged information by using a reduced chaos expansion 
with random coefficients. This expansion has recently been proposed in [13], and in the 
current context, it enables us to obtain a reduced-dimensional representation of the exchanged 
information, a key feature of which is that it can maintain segregation between sources 
of imccrtainty that stem from different subproblems. Further, we present a corresponding 
measure-transformation technique that allows stochastic expansion methods to exploit this 
reduced-dimensional representation to obtain an efficient solution of subproblems in a reduced- 
dimensional space. We show that owing to the uncertainty source segregation, requisite 
orthonormal polynomials and quadrature rules can be readily obtained by tensorization. 

The remainder of this paper is organized as follows. First, in Sec. 2, we outline the 
proposed methodology. Next, in Sec. 3, we describe the reduced chaos expansion with random 
coefficients. In Sec. 4, we provide details on the implementation. Finally, in Sees. 5 and 6, we 
demonstrate the proposed methodology by considering an illustration problem. 
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2. Proposed methodology: Dimension reduction and measure transformation 

2.1. Model problem 

This paper is devoted to the determination of the solution to a stochastic coupled model of 
the following form: 

f{u.x.C} = Q, y = h{u,^), / : R*^ X M*" X R" ^ h : W x ^ M.''" , 
9{y,v,C) = 0, x = k{v,C), g : R""" X R** X R" -S^ R^ fe : R** x R" -i^ R"" . ^' 

To avoid certain technicalities involved in infinite-dimensional representations, we assume that 
these equations are discretized representations of a stochastic model that couples two physical 
processes, two scales, two domains, or a combination of these subproblems. For instance, these 
equations may be obtained by the spatial discretization of a steady-state problem or they may 
be obtained at a single time step after the spatial and temporal discretization of an evolution 
problem. Further, we assume that the data of the first subproblem, which enter this subproblem 
as coefficients or loadings or both, depend on a finite set of uncertain real parameters denoted 
as ^1, . . . and that the data of the second subproblem depend on a finite set of uncertain 
real parameters denoted as We collect these sources of uncertainty into vectors 

^ = {£,1, ■ ■ ■ ,(,m) and ^ = (Ci, . . . , C„), which we model as random variables that are defined on 
a probability triple {Q,T,P) and which take values in R™ and R", respectively. Throughout 
this paper, we refer to ^ and C as input random variables. Lastly, we assume that the input 
random variables are statistically independent of one another, that is, we assume that ^ and ^ 
are statistically independent random variables. 

The stochastic coupled model in (1) is a general bidirectionally coupled model. The solution 
random variable u of the first subproblem depends on the solution random variable v of 
the second subproblem via the coupling random variable x. Likewise, the solution random 
variable v depends on the solution random variable u via the coupling random variable y. 

Thus, to solve this stochastic coupled model, it is necessary to find the random variables u 
and V defined on {Q,T,P) and with values in R*" and R* such that (1) is satisfied under the 
assumption that the stochastic coupled model is well posed, that is, the model admits a unique 
and stable solution. 

Finally, it should be stressed that the model problem given by (1) is identical to the one given 
in [11] and [12], except for the fact that ^ and ^ are required to be statistically independent 
here; these input random variables were allowed to be statistically dependent in [11] and [12]. 

2.2. Partitioned iterative solution 

Because a coupled model usually defines its response only in an implicit manner, the numerical 
solution to a coiiplcd model is typically obtained by using an iterative method. Here, we assume 
that iterative methods and associated solvers already exist for each subproblem. Therefore, to 
solve the coupled model, we consider a partitioned iterative method that reuses these iterative 
methods as steps in a global iterative method built around them. Let us assiime that each of the 
iterative methods is based on the reformulation of the associated subproblem as a fixed-point 
problem as follows: 

u = a{u,x,^), y = h{u,$), a : R'^ x R^o x R" ^ R'', h : R'' x R™ R''° , 

v = b{y,v,C), x = k{v,C), b:W°xWxW^W, fe : R** x R" ^ R**". ^' 
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Wc note that these equations can be obtained by setting a{u,v,^) = u — f{u,v,^) 
and b{u,v,(^) = v — g{u,v,(^), but alternative reformulations, such as those involving direct 
solutions of the subproblems or of their linear approximations, are often better adapted. We 
then consider the solution of the stochastic coupled model by a Gauss-Seidel iterative method, 
using suitable initial values u'^, v'^, and x''^ = k(v''\ C), as follows: 

This is not the only partitioned iterative method available; however, for simplicity, we employ 
only this method in this work. Although we implement the proposed methodology by using 

the Gauss-Seidel iterative method, we note that the proposed methodology can be readily used 
with other iterative methods, such as Jacobi, relaxation, and Newton methods. 

2.3. Dimension reduction 

We believe that exchanged information often resides in a considerably lower dimensional 
space than the sources of uncertainty themselves. In [11] and [12], we had thus investigated 
the effectiveness of dimension-reduction techniques for the representation of the exchanged 
information; specifically, rather than the exchanging of the coupling random variables x^ and 
and solution random variables and in their original form, we had proposed the 
approximation of these random variables by using a truncated KL decomposition as they 
pass from subproblem to subproblem and from iteration to iteration. Here, we present an 
extension of this previous work. For stochastic coupled problems whose subproblems have 
only input random variables that are statistically independent of one another, we present 
a characterization of the exchanged information by using an alternative dimension-reduction 
technique, namely, by using a reduced chaos expansion with random coefficients. This expansion 
has recently been proposed in [13], and in the current context, it can be applied to solve the 
stochastic coupled model as follows: 



(4) 



v' = b{y'-',v'-'^'^,C), x' = k{v'X\ 

qiA _ ,^£-i,d-j j^j^j ^i,e _ ({j^^i.e ;j;^-i.e) are reduced chaos expansions with 

random coefficients of = {y^,v'^~^) and = {u^^^^ ,x^^^), respectively, which read as 

The reduced chaos expansion with random coefficients is discussed in detail later. However, 
here, it is to be noted that these expansions take the form of KL decompositions with random 
basis vectors. We say that these expansions maintain segregation between the input random 
variables because each expansion captures one input random variable in the usual manner using 
the reduced random variables, while capturing the other input random variable, statistically 
independent from the first, by constructing the basis vectors as random vectors. 
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Because truncations often result in approximation errors, we use a hat superscript to 
distinguish between successive approximations determined by (3) and (4). 

We note that (5) provides a combined reduced-dimensional representation of and v^~^ in 
terms of the reduced random variables ij^ = {rj^, ...,77^), and a combined reduced-dimensional 
representation of tt^"^ and x^~^ in terms of the reduced random variables = . . , 4)- 
However, these are not the only dimension reductions that can be considered. The proposed 
methodology can be readily used with variants of this dimension reduction, such as variants 
involving separate reduced-dimensional representations of the coupling and solution random 
variables, with each representation having its own reduced random variables. 

Further, although our notations do not express the potential dependence of d and e on i, 
we note that the reduced dimensions may depend on the number of iterations. 

2.4- Measure transformation 

The successive approximations in (3) determined by the iterative method that does not involve 
dimension reduction can be constructed as random variables of the following form: 

that is, and can be constructed as transformations of | = (^i,...,Cm) and C = 
(Ci, • • • 5 Cn)- Thus, and exist in solution spaces of stochastic dimension m + n. Hence, 
implementations of (3) using stochastic expansion methods typically involve the approximation 
of the transformations that map ^ and ^ into and by finite-dimensional representations. 

Depending on whether the nonintrusive projection, embedded projection, or collocation 
method is chosen, this approximation requires the construction of orthonormal basis functions, 
quadrature rules, moment tensors, or a combination of these with respect to the joint 
probability distribution of ^ and Owing to the assumption that ^ and C are statistically 
independent, the joint probability distribution of ^ and C is the product probability 
distribution x P(^, where and P^ are the probability distributions of ^ and 
respectively. Thus, orthonormal basis functions and quadrature rules with respect to P^ x 
can be constructed by tensorization from orthonormal basis functions and quadrature rules 
constructed separately with respect to P^ and P,^; for details, the reader is referred to [5] 
and [14]. Hence, implementations of (3) typically involve the approximation of and as 

l«H-|/3|=0 
l«H-|/3|=0 

where {(fia,Oi G N"*} and {ipp,^ G N"} are orthonormal bases for the spaces of P^-square- 
integrable functions from R™ into R and of P^-square-integrablc functions from R" into R, 
respectively. These orthonormal bases are indexed by the usual multi-indices cx = (ai, . . . , am) 
in N™, with |a| = ai + . . . + a„, and /3 = (A, • ■■,Pn) in N", with |/3| = A + • • • + ^n- 
In contrast, owing to the dimension reduction, the successive approximations determined by 



</3 e R^ 



v:,^ e R^ 



(7) 
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the iterative method (4) can be constructed as random variables of the following form: 

that is, can be constructed as a transformation of ^ = (^i, . . . ,Cm) and = . . . , tg) 
and i)^ can be constructed as a transformation of 77^ = (r/f , . . . , r]^) and (J = {(i,. . . ,(n)- Thus, 
and i)^ exist in solution spaces of stochastic dimensions m + e and d + n, respectively. 
Because ^ and ^ are statistically independent and because the reduced chaos expansions 
with random coefficients given by (5) imply that is a transformation of only ^ and r/^ is a 
transformation of only ^, the random variables ^ and and, likewise, the random variables r)^ 
and C are statistically independent of one another. Thus, orthonormal basis frmctions and 
quadrature rules with respect to the joint probability distribution x of ^ and and 
those with respect to the joint probability distribution x of rf^ and C can be obtained by 
tensorization from orthonormal basis functions and quadrature rules constructed separately 
with respect to the probability distributions P^ and P/ and from those constructed separately 
with respect to the probability distributions P,^ and P^, respectively. Hence, implementations 
of (4) can exploit the dimension reduction to approximate and as follows: 

|a|+|c|=0 
|7| + I/3|=0 

where {T^,/3 g N^} and {r^,7 € N''} are orthonormal bases for the spaces of P/-square- 
integrable functions from M*^ into R and of P^-square-integrable functions from M*^ into R, 
respectively, indexed by the usual multi-indices <; = {<;i, ...,<;(,) in N^, with |<j| = Ci + • . . + <?e> 
and 7 = (71, . . . , 7d) in N'', with \^\ = + . . . + 7^. 

We say that this approach involves a measure transformation because it follows from (5) that 
the probability distributions of the reduced random variables and t]^ are transformations of 
the probability distributions of the input random variables ^ and as described later. 

Throughout this work, we employ orthonormal bases that consist of polynomials and we 
thus refer to them as polynomial chaos bases. 

We note that although our notations do not express the potential dependence of p and q on 
the subproblem or i, the subsets of the basis functions used to construct the finite- dimensional 
representations may depend on the subproblem and the number of iterations. 

2.5. Effectiveness of the reduced chaos expansions with random coefficients 

In the reduced chaos expansions with random coefficients, the statistical independence of ^ 
and and that of 77^ and which are a consequence of the statistical independence of ^ and ^, 
are significant because these statistical independences enable polynomial chaos and quadrature 
rules with respect to the joint probability distribution P^ x P^ and those with respect to the 
joint probability distribution P^ x P^ to be obtained by tensorization from polynomial chaos 
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and quadrature rules obtained separately with respect to and and from those obtained 
separately with respect to P,^ and P^, respectively. 

In contrast, even though the methodology presented in [11] and [12] applies to more general 
stochastic coupled problems whose subproblcms need not have only input random variables 
that are statistically independent of one another, the use of KL decompositions for the 
representation of the exchanged information does not allow the statistical dependence of ^ 
and C to be exploited for the construction of the requisite polynomial chaos and quadrature 
rules by using tensorization. In general, the use of KL decompositions leads to reduced- 
dimensional representations of the exchanged information, whose reduced random variables — 
let us still denote them as and r)^ — are each a transformation of both input random 
variables ^ and (J, and whose basis vectors are deterministic vectors. As discussed in [11] 
and [12], the implementation of the associated measure transformation leads, in this case, to 
the construction of polynomial chaos and quadrature rules with respect to the joint probability 
distributions P^^ (of ^ and l^) and P^^ (of t]^ and C) because the random variables ^ and 
and the random variables rj^ and ^ are then statistically dependent on one another. 

In many applications, the components of the input random variables ^ and Q are statistically 
independent "labeled" random variables; then, polynomial chaos and quadrature rules with 
respect to P^ and P(^ can be readily constructed by tensorization of univariate polynomial chaos 
and quadrature rules read from tables in the literature [4, 5, 14-18]. Well-known examples 
include the use of Hermite and Legendre polynomial chaos and of Gauss-Hermite and Gauss- 
Legendre quadrature rules with respect to Gaussian and uniform probability distributions. 

However, the reduced random variables and rj^ associated with reduced chaos expansions 
with random coefficients or KL decompositions are usually statistically dependent and 
not "labeled." Hence, usually, polynomial chaos and quadrature rules with respect to the 
probability distributions P/ and P^ and, likewise, polynomial chaos and quadrature rules with 
respect to the joint probability distributions P^^^ and P^^^ cannot be read from tables in 
the literature; they should be computationally constructed. 

Now, the computational cost associated with the computation of polynomial chaos and 
quadrature rules with respect to a given probability distribution can be expected to increase 
with an increase in the dimension of the space on which this probability distribution is defined. 
Thus, for stochastic coupled problems whose subproblems have only input random variables 
that are statistically independent of one another, the use of reduced chaos expansions with 
random coefficients for the representation of the exchanged information has the advantage 
of reducing the computational effort associated with the implementation of the measure 
transformation since only the computation of polynomial chaos and quadrature rules with 
respect to the probability distributions of the reduced random variables is necessary. 

2. 6. Effectiveness of the proposed methodology 

The main feature of the proposed methodology is that it provides a solution of the subproblems 
in a reduced-dimensional space when the exchanged information has a low effective stochastic 
dimension; specifically, a solution in a reduced-dimensional space is obtained when the reduced 
dimensions can be selected such that d < m and e < n, while sufficient accuracy is maintained 
(refer to (6) and (8)). This is a significant benefit since stochastic expansion methods suffer 
from a curse of dimensionality: their computational cost increases quickly with an increase in 
the stochastic dimension. The proposed methodology breaks the curse of dimensionality by 
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inhibiting the increase in stochastic dimension when information is exchanged. 

Finally, we note that the proposed methodology can be adapted to meet various requirements 
of specific applications. For instance, one could only use a dimension-reduction technique to 
represent exchanged random variables that arc of low effective stochastic dimension or one 
could only implement a measure transformation to solve subproblems whose computational 
cost would thereupon be lowered. 



3. Reduced chaos expansion with random coefficients 

5.1. Problem setting 

Let ^, and q be three random variables defined on the probability triple (0, T, P) with values 
in the Euclidean spaces R™, M", and , respectively; let their probability distributions be P^, 
P^, and Pq. Further, let the random variables ^ and C be statistically independent; hence, the 
joint probability distribution of ^ and Q is the product probability distribution x P^. Lastly, 
let the random variable q be a transformation of ^ and ^ given by 

q{e)^q{m,m) (10) 

and let it be of the second order: 

/ / \\q{i,Q\?dP^dP^<+^, (11) 

where ||-|| denotes the Euclidean norm. Under these assumptions, we construct below a reduced 
chaos expansion with random coefficients of q. 

3.2. Interpretation 

In the context of the proposed methodology, we consider q as a random variable that passes 
between subproblems or iterations during the determination of the solution to a stochastic 
coupled model. With reference to (5), at a specific iteration, q could collect the components 
of and v^~^ with w = ro + s ov those of u^~^ and x^~^ with w = r + sq. 

3.3. Chaos expansion with random coefficients 

Let a € N"*} and {'4>i3,P € N"} be polynomial chaos bases for the spaces of P^-square- 
integrable functions from R™ into R and of P^-square-integrable functions from R" into R, 
respectively, indexed by the usual multi-indices a = (ai, . . . , a^) in and /3 = (/3i, . . . , /?„) 
in N". Throughout this paper, we assume that ipQ = l and Vo = Ij hence, we have 

/ V<.{£.)dPi=G. a^O, and / ^p{OdPi^=Q, f3 ^ 0. (12) 

Owing to the statistical independence of ^ and ^ and, thus, to the fact that their joint 
probability distribution is the product distribution P^ x P(^, the random variable q has the 
following chaos expansion related to the tensor product of the aforementioned spaces: 

9(^'C)= E E 9«/3 ^0(0^-/3(0, / QiiX)'PUi)MC)dP^dPc; (13) 
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9(^,0 



E 



E iisiOMO, 



(14) 



here, both sets of sources of uneertainty, that is, ^ and ^, are captured in the usual manner 
using polynomial chaos, and the coefficients q^/j are constructed as deterministic vectors in M}" . 

Following the definitions given in [13], the random variable q has the following two 
corresponding chaos expansions with random coefficients: 

9c. (C) = E 9«/3V'/3(C), 

/3GN" 

here, each chaos expansion with random coefficients captures one set of sources of uncertainty 
in the usual manner using polynomial chaos and the other set by constructing the coefficients as 
random variables. It should be emphasized that because ^ and ^ are statistically independent 
and because the expansions given by (14) imply that the random coefficients are 
transformations of only ^ and the random coefficients are transformations of only 4, 
the random coefficients q„ are statistically independent of ^ and, likewise, the random 
coefficients are statistically independent of ^. 

Below, we reduce the second of the two chaos expansions with random coefficients in (14) 
to obtain the corresponding reduced chaos expansion with random coefficients. We note that 
the first of the two chaos expansions with random coefficients in (14) can be readily reduced 
in a similar manner; however, for the sake of brevity, this reduction is not shown. 



3.4- Reduced chaos expansion with random coefficients 

Let q be approximated by a chaos expansion with random coefficients truncated at total 
degree p as follows: 



E 10(0 MO, 

1/31=0 



(15) 



wherein the random coefficients are defined by (14). Then, according to [13], an associated 

reduced chaos expansion with random coefficients is obtained by approximating the random 
coefficients in this expansion by using a truncated KL decomposition. 

Let n = card({'!/'/3, |/3| < p}) = (n + p)\/n\/p\ denote the number of terms in the chaos 
expansion with random coefficients given by (15). Then, let the random coefficients {q^, |/3| < 
p} of qP be collected in a random variable with values in 



(16) 



where /3i, . . . are the multi-indices of the random coefficients contained in {q^, |/3| < p} 

arranged in a sequence. With reference to (12), the mean vector of this random variable is 
the {fj, X t(;)-dimensional vector that collects the blocks q/s = qp{^)dP^, and it can be 
expressed as a function of the coefficients as follows: 







Sa6N 



















.9/3^. 







(17) 
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Likewise, the covariancc matrix is the (p x ?j7)-dimensional square matrix that collects the 

blocks Cq^q^ = (9/3(0 — Qfi) ~ Q^) dP^, and it can be expressed as a function of 

the coefficients q^p as follows: 



C„ 



(18) 



Let W be a tw-dimensional, square, symmetric, and positive definite matrix, termed the 
weighting matrix. Because the abovementioned covariance matrix is symmetric and positive 
semidefinitc, the solution of the generalized cigcnproblcm 

rT, 



W^Cq^..,.W ... W^Cq W 









W 








A. 




W 







(19) 



provides jjLXw eigenvalues, Ai > A2 > . . . > A^xtu > 0, and i^xw eigenvectors, (f)^,. . . 
These eigenvectors form a W-weighted orthonormal basis of IR''^'", that is. 



tJiXW 



w 






w 




A. 



(20) 



where 5ij is the Kronecker delta; Stj is equal to 1 if i = j, and it is otherwise. Then, the KL 
decomposition of the vector given by (16) is obtained as 



"9/3,(0" 






flXW 



















where the reduced random variables r]j are random variables with values in M such that 
1 



(21) 



(9/3,(0-9/3,)'^ ••• (9/3 (0 -9/3 J ^ 



w 




'^k 


w 




A. 



These reduced random variables r]j are zero-mean and uncorrelated: 

By substituting (21) in (15), the following representation of 9^ is obtained: 

flXW 

9^(0 c) = 9^(0 + E V^^^ ii)<f>'''{c)- 



(22) 

(23) 

(24) 

(25) 
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This representation takes the form of a KL decomposition that captures one input random 
variable, ^, in the usual manner using the reduced random variables and the other input 
random variable, C by constructing the basis vectors as random vectors, that is, 

r(C) = E ^opMci (26) 

|/3|=0 

0^>(C) = (27) 

|/3|=0 

By truncating (25) after d terms, we obtain a reduced chaos expansion with random 
coefficients q^''^ of q as follows: 

d 

q^'\i, C) = 9^(0 + E ^/^imiiW'^iO- (28) 
Because of the orthonormality properties (20) and (24), the truncation error satisfies 

q^{iX)-q'''\iX)fwdP^dPc= E ^i' (29) 

where ||-||^ is the VT-weighted norm on W" such that = for any p in W" . 

Thus, the accuracy of the reduced chaos expansion with random coefficients in (28) can be 
improved systematically by increasing the number of terms that are retained. 

3.5. Optimality 

Let (8) M'" be the space of random variables with values in that are representable as 
a chaos expansion in ^ of maximum total degree p. Then, for any basis {e-'^'^, . . . , e''^"''^} 
of pP (Si with e^'P = J2^fi\=o ^^V'/S that is VF-weighted orthonormal in that 

p 

|/3|=0 

the random variable qP can be expanded as follows: 

9^(^,0= g^(C)+E E {Qi3{^)-Ql3fWe^]e^'P{0. (31) 

j = l ^|/3|=0 ^ 

The approximation error = J2'j=i Sf/3|=o(9/3 ~ qi3)'^We-'pe^ P introduced because of the 
truncation of this expansion after d terms can be readily shown to satisfy 

\\e,{^,C)\\wdPidPc^ E E E i^pfw^C.^i.Wel. (32) 

j=d+l |/3|=0 1^1=0 
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Expression (32) indicates that the redueed ehaos expansion with random coefheients is 
optimal in that from among all expansions with d random basis vectors that are VF-weighted 
orthonormal, it minimizes the W-weighted mean-square norm of the approximation error: 

Expression (33) indicates that the reduced chaos expansion with random coefficients has 
an optimality property that is similar to that of the KL decomposition; for details on the 
optimality property of the KL decomposition, the reader is referred to [4] . However, we note 
that whereas the KL decomposition achieves optimality among all expansions of a given finite 
length that capture all sources of uncertainty using the reduced random variables and that 
construct the basis vectors as deterministic vectors, the reduced chaos expansion with random 
coefficients achieves optimality among all expansions of a given finite length that capture 
one set of sources of uncertainty using the reduced random variables and the other set by 
constructing the basis vectors as random vectors. 

3.6. Concluding remarks 

The reduced chaos expansion with random coefficients described here can be viewed as an 
adaptation of the one introduced in [13]; the construction presented here features a weighting 
matrix, whereas the one given in [13] does not. This weighting matrix is particularly useful 
for the construction of a reduced-dimensional representation of a random variable that solves 
a space-time discrctized stochastic model. Then, by choosing the weighting matrix as the 
Gram matrix of the discretization basis, we obtain a reduced-dimensional representation that 
is consistent with Hilbertian projections in the Hilbert space of random variables to which that 
random variable belonged prior to the discretization of the stochastic model. Therefore, the 
representation satisfies the optimality condition (33) in a weighted norm that is consistent with 
the norm that was relevant to that random variable prior to the discretization of the stochastic 
model; for more details, the reader is referred to [11]. Nevertheless, the reduced chaos expansion 
with random coefficients introduced in [13] is recovered by setting the weighting matrix equal 
to the identity matrix. 



4. Implementation 
4.1. Reduced chaos expansion with random coefficients 

In this section, we describe the implementation of the reduced chaos expansion with random 
coefficients of a random variable that is represented by a chaos expansion involving only 
polynomial chaos whose total degree does not exceed a given value. We show how this 
implementation in turn naturally provides a representation of the reduced random variables 
by a chaos expansion involving only polynomial chaos whose total degree does not exceed the 
given value. Specifically, we adopt the notations used in Sees. 2 and 3, and we construct a 
reduced chaos expansion with random coefficients of a second-order random variable with 



mm 

{ei''',...,e'''»'}epP(8)M' 



/ / hifwdP^dPc- (33) 
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values in that is represented by the following chaos expansion: 

g^(C,C)= E (0^-/3(0, g«/3eR-. (34) 

l«|+|/3|=0 

Then, can be expressed by the following chaos expansion with random coefficients: 

P P-|/3| 

9^(^,0 = E 9r""(^)^/3(C), Qr""(^) = E 9«/3^a(^). (35) 

|/3|=0 |a|=0 

Owing to the truncation of the chaos expansion of given by (34) at total degree p, we obtain 
a representation of each random coefficient g^ ''^' that has the form of a chaos expansion 
involving only polynomial chaos up to total degree p — |/3|. Because of the orthonormality 
of the polynomial chaos (pa, the mean vectors and cross-covariance matrices of the random 
coefficients are immediately obtained as follows: 

p-max(|/3|,|^|) 

Q/3 = Qo/3 and Cg^g^ = E ^"/s^a^' (^6) 

l«|=i 

where owing to the truncation of the chaos expansion of q^ given by (34) at total degree p, 
cross-covariance matrices Cg^g^ for which |/3| > p — 1, or |/3| > p — 1, or both, vanish. Then, 

the solution of the generalized eigenproblem ^^^i^q^^*^^ ~ ^j^^js' 1/^1 ^ P~ 1) 

provides the eigenvalues and associated cigenmodcs required to construct a reduced chaos 
expansion with random coefficients qP''^ of q^ as follows: 

d 

Q^'^it = g^(C) + E VaJtjJ (^)0^>-^(C), (37) 

where owing to the truncation of the chaos expansion of q^ in (34) at total degree p, 
the basis vectors qP and 0^'^"^ are represented by chaos expansions q^ ~ S/3=o^/3V'/3 
and <p^'P~^ = J2]^/3\=o 't^isi'P involving polynomial chaos up to total degrees p and p — 1. The 
reduced random variables rjj are random variables with values in M such that 

^(0 = ^ E (9r""(^) -^/3)^W^0^, (38) 

and they are zero- mean and uncorrelated. By substituting (35) in (38), we obtain the 
representation of each reduced random variable as a chaos expansion: 

p ^ p-i 

= E ^i.°^"(C) with r?,-,„ = ^ E llfiW'P^is- (39) 

|a|=l V^i l/3|=o 

Expression (39) indicates that the reduced chaos expansion with random coefficients of a 
random variable that is represented by a chaos expansion involving only polynomial chaos 
whose total degree does not exceed a given value naturally provides a complete probabilistic 
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characterization of the reduced random variables by a chaos expansion involving only 
polynomial chaos whose total degree does not exceed the given value. 

We note that while this section provides details on the implementation of the reduced chaos 
expansion of a random variable that is represented by a chaos expansion that is truncated at a 
given total degree, the proposed implementation can be readily extended to random variables 
represented by chaos expansions involving other subsets of polynomial chaos, such as those 
obtained by full tcnsorization; for examples, the reader is referred to [19]. 

Finally, we note that the methodology proposed in Sec. 2 for obtaining the solution 
to stochastic coupled problems provides a representation of the solution and coupling 
variables associated with the subproblcms as chaos expansions in terms of combinations of 
the input random variables ^ and ^ and the reduced random variables and r]^ of the 
reduced-dimensional representations of the exchanged information (refer to (8)). Thus, when 
these reduced random variables arc represented by chaos expansions in terms of the input 
random variables ^ and C themselves (refer to (39)), the proposed methodology requires the 
construction of reduced chaos expansions with random coefficients of random variables that are 
represented by compositions of chaos expansions. The implementation of such reduced chaos 
expansions with random coefficients falls within the scope of the implementation mentioned 
above because the composition of two chaos expansions can always be written equivalently as 
a chaos expansion of the form of (34), albeit as one that is truncated at a higher total degree. 

4-2. Construction of the requisite polynomial chaos 

Depending on the stochastic expansion method that is chosen, the approximation of the 
solutions to the subproblcms by chaos expansions in terms of combinations of the input 
random variables ^ and C and the reduced random variables and r/^, proposed in Sec. 2, 
requires the construction of polynomial chaos (discussed in this section) and the construction of 
moment tensors, quadrature rules, or both (discussed in the next section). Here, we construct 
polynomial chaos with respect to the product probability distribution x P^, which is 
defined on the Euclidean space and which involves the probability distribution of 

the reduced random variable rj^ with values in R'^ and the probability distribution P(^ of the 
input random variable ^ with values in R". Clearly, polynomial chaos with respect to P^ x P/ 
can be constructed in a similar manner; however, this construction is omitted here for brevity. 

4-2.1. Polynomial chaos with respect to P^. In many applications, the components of the 
input random variables are statistically independent "labeled" random variables. Thus, let 
us assume that the probability distribution of the input random variable C is the product 
probability distribution P^ = Pq-^ x . . . x Pq^ consisting of "labeled" univariate probability 
distributions P.^^ , . . . , P.^^ . Then, for each dimension j = l,...,n, a basis G I^} of 

polynomial chaos of increasing degree for the space of P^^. -square-integrable functions from R 
into R can be readily read from tables in the literature, and a corresponding basis {V'/a, /3 G N"} 
of polynomial chaos of increasing total degree for the space of P|j-square-integrable functions 
from R" into R, indexed by the multi-indices (3 in N", can be readily obtained by tcnsorization 
by setting t/'/3(C) = V'^i(Ci) x • • • x V'^^(Cn); for details, the reader is referred to [4, 5, 14, 17]. 

4.2.2. Polynomial chaos with respect to P^. The reduced random variables associated with 
a reduced chaos expansion with random coefiicients are usually statistically dependent and 
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not "labeled." Hence, polynomial chaos with respect to the probability distribution of the 
reduced random variable rj^ usually cannot be read from tables in the literature, and they 
should be computationally constructed. In this paper, we use the procedure presented in [12], 
wherein wc first arrange the multivariate monomials rj'^ = rfl^ x . . . x ryj'' in a sequence of 
increasing total degree and then orthonormalize this sequence by means of the Gram-Schmidt 
method; we use the inner product associated with to obtain the requisite polynomial chaos. 
Let {r^,7 G N''} denote the polynomial chaos basis thus obtained for the space of P^-square- 
integrable functions from into R, indexed by the multi-indices 7 in N''. 

4-2.3. Polynomial chaos with respect to P^ x P^. Once the polynomial chaos bases {r^,7 G 
N'^} and {^jp^j3 G N"} for the spaces of P^-square-integrable functions from into M and 
of P^-square-integrable functions from R" into M, respectively, are available, a corresponding 
polynomial chaos basis for the space of P^ x P^ -square- integr able functions from into K is 
readily obtained by tensorization as {r^?/'/3j 7 G '^'^^ /3 € N"}, as discussed in Sec. 4.2.1. Hence, 
any (P^ x P^)-square-integrablc function / from M''+" into R can be expanded as follows: 



/('7^C)= E E /-T/3r^(^')V'/3(C), A/3=/ / f{n'xK{'n')MC)dP^dPc. (40) 



Then, various finite subsets of polynomial chaos can be selected to construct chaos expansions 
of finite length; for examples, the reader is referred to [19]. In this work, we use finite subsets 
of the form {r^^^g, I7I + < q\ to construct chaos expansions of the form of (9), involving 
only polynomials whose total degree does not exceed a given value q. 

4.3. Construction of the requisite quadrature rules 

In this section, we consider the construction of a family of quadrature rules for integration 
with respect to the product probability distribution P^ x P^, which is defined on the Euclidean 
space R''+"' and which consists of the probability distribution P^ of the reduced random 
variable r]^ with values in R'' and the probability distribution P^ of the input random variable C, 
with values in R". Clearly, a family of quadrature rules with respect to P^ x P/ can be 
constructed in a similar manner; however, this construction is omitted here for brevity. 

4.3.1. Quadrature rules with respect to P^. As in the previous section, let us assume that P^ 
is the product probability distribution given by P^ = Pq^ x . . . x P^^ and consisting of the 
"labeled" univariate probability distributions Pq^ , . . . , Pq^ . Then, for each dimension j = 
1, . . . , n, families of quadrature rules of increasing accuracy for integration with respect to Pq. 
can be readily constructed and often read from tables in the literature. While quadrature 
rules for integration with respect to univariate probability distributions can be obtained 
following various approaches [15, 16, 18], we have used, in this study, Gaussian quadrature 
rules for integration with respect to univariate probability distributions. Thus, for each 
dimension j = 1, . . . , n, let Q^. be the level-A Gaussian quadrature rule that allows the integral 
of any continuous function / from R into R with respect to P^^ to be approximated by a 
weighted sum of integrand evaluations as follows: 



-ceN-* /Sew 




A 




(41) 
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It is well known that a Icvcl-A Gaussian quadrature rule allows integrals of univariate 
polynomials up to degree 2A — 1 to be calculated exactly. 

Because the components of ^ = (Ci, ■ • ■ ,Cn) are statistically independent, quadrature rules 
for integration with respect to the multivariate probability distribution can be readily 
synthesized from the aforementioned quadrature rules for integration with respect to the 
univariate probability distributions P^^ , . . . , Pf„ by tensorization. While various approaches are 
available to obtain such tensorized quadrature rules [18], in this work, we use fully tensorized 
quadrature rules if the stochastic dimension is low, say n < 5, and sparse-grid quadrature rules 
otherwise. A level-A fully tensorized quadrature rule that allows the integral of any continuous 
function / from R" into M with respect to P^ to be approximated by a weighted sum of 
integrand evaluations is obtained as follows: 

A A 

(Q^, ® . . . ® Ql){f) = E ■ ■ ■ E • • X • • • X <K- (42) 

fel = l hn=l 

Further, a level-A sparse-grid quadrature rule is obtained as follows: 

X<\l\<\+n-l ^ ^ / I 1/ 

Where qJ;^ ® . . . ® qJ- )(/) = ■ ■ ■ ^t=i fiCiX'- ■ ■ ' Cl)^i':^, x . . . x z;^^. Various 

approaches are available to select A^^. as a function of Ij in (43), for example, the classical 
Smolyak approach with A;^ = 2'^ — 1 and the so-called slowly increasing approach with A;^ = Ij. 
In this work, we use the Smolyak approach. 

By the argument given in [20], it can be shown that if each univariate quadrature rule of 
level A;^ allows integrals of univariate polynomials of at least up to a degree of 21 j — 1 to be 
calculated exactly, then the sparse-grid quadrature rule of level A given by (43) allows integrals 
of multivariate polynomials of at least up to a total degree of 2A — 1 to be calculated exactly. 

4-3.2. Quadrature rules with respect to P^. The reduced random variables associated with 
a reduced chaos expansion with random coefficients are usually statistically dependent and 
not "labeled." Hence, usually, quadrature rules for integration with respect to the probability 
distribution P^ of the reduced random variable r}^ cannot be read from tables in the literature, 
and they should be computationally constructed. 

In this work, we use the procedure proposed and described in detail in [12]. This procedure 
is applicable, provided that the reduced random variable is itself characterized by a chaos 
expansion in another random variable whose probability distribution is such that a quadrature 
rule for integration with respect to this probability distribution can be easily constructed, for 
example, r]^ = r)^'^{^) = X]|'a|=i vLfai^) (refer to (39)). In this procedure, first, a quadrature 
rule for integration with respect to P^ is constructed. Next, a change of variables is carried 
out to obtain a quadrature rule Q^{- orj^'P) for integration with respect to P^. Finally, a subset 
selection algorithm is used to obtain an embedded quadrature rule, 

Qii\f) = Y:f{v',')wj:\ (44) 
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that uses a subset of the nodes of the quadrature rule Q^{- o tj^^p) as nodes and that has level A 
in that it allows integrals of multivariate polynomials of at least up to a total degree of 2A — 1 
to be calculated exactly. 

4-3.3. Quadrature rules with respect to x P^. Once A-indexed families of quadrature 
rules Qlj^ and for integration with respect to P^ and P,^ are available, corresponding 
quadrature rules for integration with respect to P^ x P^ can be synthesized by tensorization. 
As mentioned in Sec. 4.3.1, various approaches are available to obtain tensorized quadrature 
rules. In this work, we use the sparse-grid construction. A level-A sparse-grid quadrature rule 
that allows the integral of any continuous function / from R''+" into M with respect to P^ x 
to be approximated by a weighted sum of integrand evaluations is obtained as follows: 

QfL)(f) = E (-l)^+^-('=+')(Q^'''= ® Qc')(f)^ (45) 

A<fe-l'i<A+l 

where (Qij^- ^ Q'^'){.f) = Ef^^Efi fi^l^" ^ ^i"' {(^f "-'4'^')' ^ ^ ^ ^ 

vj^^} and {(C;-' ) ')) 1 ^ ^ ^ ^A;} are the nodes and weights of Q^j^'' and Q^', respectively. 
As mentioned in Sec. 4.3.1, various approaches are available to choose A^ and Aj as a function 
of k and I. In this work, we use here the slowly increasing approach with Afc = fc and Xi = I. 

Using the argument given in [20] , it can be shown that if the component quadrature rules 
of levels Afc and A; allow integrals of multivariate polynomials of at least up to total degrees 
of 2Afc — 1 and 2 A; — 1, respectively, to be calculated exactly, then the sparse-grid quadrature 
rule of level A given by (45) allows integrals of multivariate polynomials of at least up to a 
total degree of 2A — 1 to be calculated exactly. 



5. Realization for a multiphysics problem 

5.1. Problem formulation 

[Figure 1 about here.] 

We consider the stationary transport of neutrons in a one-dimensional reactor with 
temperature feedback [21]. Let the reactor occupy an open interval ]0, L[ (Fig. 1). The problem 
involves finding the temperature T and neutron flux $ such that 

(46) 

under homogeneous Neumann boundary conditions. The first term on the left-hand side of the 
heat subproblem represents heat conduction, and the second term represents the transmission 
of heat to the surroundings; further, the right-hand side represents a distributed heat source 
proportional to the neutron flux. The flrst term on the left-hand side of the neutronics 
subproblem represents neutron diffusion, and the second term represents the net effect of the 
absorption and generation of neutrons; further, the right-hand side represents a distributed 
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neutron source. The coefficients k and h arc the licat conductivity and licat transmittance, 
respectively; the temperature T^o is the ambient temperature; and u and E{ are the number 
of neutrons and the energy released per fission reaction, respectively. The coefficients D, E^, 

and Ef arc the neutron diffusion constant, fission cross section, and absorption cross section, 
respectively; the dependence of these coefficients on the reactor temperature is as follows: 



D{T{X)) = Aef^ Mn^)) = ^a,refy Sf (T(a;)) = E,,,,^ (47) 

5.2. Deterministic weak formulation 

Let H = H^{]0,L[) be the space of functions that are sufficiently regular to describe the 
solutions to the heat and neutronics subproblems. Then, the weak formulation involves 
finding T and $ in iJ such that 



/ k- — -dx+ / h(T-T^)Sdx= / EfJ:i(T)^Sdx, G H, 

Jo dx dx Jo Jo 

/ D(T)- — -dx+ / (Ea(T)-z/Ef(T))$*da;= / s^dx, e H. 

Jo dx dx Jo \ J Jo 



(48) 



5.3. Random thermal transmittance and absorption cross section 

We incorporate uncertainties by modeling the thermal transmittance and absorption cross 
section as random fields {h{x, < x < L} and {Ea,ref(a;, Q),^ < x < L} such that 



Sa,ref(a;, C) = ^a.ref (l + ^Cj 4^ {x)) , 

V 7 = 1 ^ 



(49) 



(50) 

where the random variables and Q are statistically independent uniform random variables 
defined on the probability triple (0,T, P) and with values in [—1,1]; thus, the random 
variables >/3Ci ^^'^ V^Cj s-r^ uniform random variables with unit standard deviation. Further, 
Xj and ^ {\j and (^ ) are the eigenvalues and eigenmodes of the eigenproblem C (</)■' ) = Xjcj)^ 
{C{4>^) = Xj<t>')\ C and C are covariance integral operators with the following kernels: 

Here, the parameters ah and as are the spatial correlation lengths of the thermal transmittance 
and absorption cross section random fields, respectively. Clearly, the random fields thus 
obtained arc such that the random variables h(x,^) and Ea.rof(2;,C) have mean values h 
and Ea.ref and coefficients of variation Sh and fc, respectively, at every position x, at least 
when the approximation errors introduced because of the truncation of the expansions after m 
and n terms, respectively, are not taken into account. 
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5.4- Stochastic weak formulation 

The weak formulation of the stochastic problem involves finding random variables T and $ 
on {e,T,P) with values in if such that with ft (^) = h{-,^) and Sa(r,C) = Sa,ref(-, C) V^ref/r, 

/ ''1-1-^^+ / h{$){T - Too)Sdx = / EiT,{{T)^Sdx, G H, 

Jo d.^ d,x Jq Jq , , 

/ DiT)- — -dx+ (Ea(T,C)-!^Sf(T))$*dx= / s'i/dx, V* e il. 
io dx dx Jq \ J Jo 

5.5. Discretization of space 

The finite element (FE) method is used for the discretization of space. The domain [0, L] is 
meshed using r — 1 elements of equal length. Lc;t Ni, . . . , be a basis of element-wise linear 
shape functions such that Nj is equal to 1 at the j-th node and at all other nodes. Using this 
basis, the random temperature T and random neutron flux $ are approximated as follows: 

r 

T-{x) = J2TjNj{x), TjGR, 

T (54) 
^^{x)=Y.^iN^{x), $,eM. 
i=i 

Then, the FE discretization of the stochastic weak formulation in (53) involves the 
determination of the random vectors T = {Ti,...,Tr) and $ = ($i,...,$r) defined 
on {Q,T,P), with values in W, which collect the nodal values of the random temperature 
and random neutron flux such that 

[K + if($)]r = g(*,T), 

[D{T) + M{TX)]^ = s. ^ ' 

Here, K, H{^), D{T), and M{TX) are r-dimensional matrices, and q(#,T) and s are r- 
dimensional vectors such that 

SjKS2= r k^^dx, (56) 
Jo dx dx 

SjH{^)S2= r h{^)S{S^2dx, (57) 

*-i>(r)*.^^^(r')^ffi^.., (58, 

^JM{T,C)^2 = (Sa(T^C) -z^Sf(T'-))*^*^dx, (59) 



S'^q{T,^)= / E{T,i{T'')'^''S''dx+ / hT^S^'dx 
Jo Jo 



(60) 



S'^s= [ s^^'dx. (61) 
Jo 
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5.6. Reformulation as a realization of the model problem 

The aforementioned illustration problem can be reformulated as a particular realization of the 
general model problem introduced in Sec. 2: 

r = a(T,*,^), a:M''xM''xM'"^M^ 

* = 6(T, C), 6 : M'^ X M" W, ^^'^^ 

where a(T,*,0 = [K + H{i)]-^q{i^,T) and 6(T,C) = [D{T) + M{TX)]-^s. This 
reformulation indicates that the illustration problem is a simplified realization of the model 
problem. This is inferred from two features. First, the neutronics subproblem admits a direct 
solution that does not require iteration. Second, the neutronics and heat subproblems are 
coupled directly through their solution variables rather than intermediate coupling variables. 

5.7. Dimension reduction 

Now, we will demonstrate the proposed methodology by approximating the random 

temperature by a reduced chaos expansion with random coefficients as it is communicated from 
the heat subproblem to the neutronics subproblem. At iteration let the random temperature 
be represented by the following chaos expansion: 

T'^'{iX)= E T'c.i3^o.{i)M0^ Kp^^'- (63) 

|c.|+|/3|=0 

Then, the random temperature has the following chaos expansion with random coefficients: 

T'^^iiX) = E ^^'"""(OV'MC), = E (64) 

|/3|=0 |c.|=0 

This expansion involves the representation of each random coefficient as a chaos expansion 
involving only polynomial chaos up to total degree p — |/3|. The mean vectors and cross- 
covariance matrices of the random coefficients are then given by 

p-max(|/3|,|/3|) 

%=fi, and Cl^^^= E Ti,{fi-^)\ (65) 

\a\=l 

where cross-covariance matrices Cq^q^ for which |/3| > p — 1, or |/3| > p — 1, or both, vanish. 
Let the r-dimensional square matrix W be the Gram matrix of the FE basis, that is. 



W 



{Nr,N,)H ... {Nr,Nr)H 



(66) 



where the inner product (•, ■) h is such that (^i, 5*2)^ = J^f SiS2dx + f^ [dSi / dx){dS2 / dx)dx 
for any pair of functions 5*1 and 5*2 in H. Then, the solution of the generalized 
eigenproblem X]fai-n ^^^^^ ^ VF^p^^ = X^Wcj)-'/, |/3| < p — 1, provides the eigenvalues 
and eigenmodes required to construct a reduced chaos expansion with random coefficients. 



T^,P,d^^^ C) = T^'^(C) + E \/^^ ^'fi^) '/''■'■""'(C). (67) 



^ V 1 '1 
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where the basis vectors T^'^ and 0-''^'^ ^ are represented by the chaos expansions T^'P = 
'^^^^Tpip/j and (f)-''^'^^^ — X]f/3|Lo '^^^^/S involving polynomial chaos up to total degrees p 
and p—1- The reduced random variables ry^'^ are random variables with values in R, such that 

rjf'i^) = 47 £ (T^^'^'^i^) - nfWcfP^', (68) 

and they are zero-mean and uncorrelated. By substituting (64) in (68), we obtain the 
representation of each reduced random variable as a chaos expansion: 

= E vUvc.{i) with = 4= E iTL,rWci>f, (69) 

l«|=i 

thus completely characterizing the reduced random variables as a chaos expansion. 

We note that the random neutron flux, in principle, could also be represented by a reduced 
chaos expansion with random coefficients as it passes from the neutronics subproblem to the 
heat subproblem. However, this extension is omitted here. 

The present implementation can be expected to be well adapted to problems wherein the 
stochastic dimension of the random thermal transmittance field is moderate or high, say m > 5, 
but the stochastic dimension of the random absorption cross section field is low, say n < 5. 
Then, the reduction of the random temperature has the potential to lower the number of 
sources of uncertainty that enter the neutronics subproblem and to thus provide a more efficient 
solution of the neutronics subproblem in a reduced-dimensional space. In contrast, a reduction 
of the random neutron flux would not have the potential to significantly lower the number of 
sources of uncertainty that enter the heat subproblem; therefore, it is not implemented. 

5.8. Measure transformation 

While ^ and C are necessarily the sources of uncertainty that enter the heat subproblem, the 
aforementioned representation of the random temperature by a reduced chaos expansion with 
random coefficients allows rj^'P and to be construed as the sources of uncertainty that enter 
the neutronics subproblem. Then, the proposed methodology leads to the approximation of 
the random temperature and neutron flux by chaos expansions as follows: 

r''^(^,C)= E Ti^MOMC), 

|a| + |/3|=0 

*^'«(,7''^C)= E *^^r^('7''^)V'/3(C), 

|7l + l/3|=0 

that is, we approximate the random temperature by a chaos expansion in ^ and ^ and the 
random neutron flux by a chaos expansion in t]^'P and C- 

We select {(fa, ct G N™} and {ippj/S G N™} as normalized Lcgendre polynomials. Following 
the approach given in Sec. 4, we construct the polynomial chaos {r^,7 e N''} at each 
iteration using the method given in [12]. Further, we select the quadrature rule for integration 
with respect to the joint probability distribution of ^ and ^, whose nodes and weights we 



Tip e M^ 



(70) 
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denote as {{^k^ Ck^'^'k), 1 < < N}, to be a sparse-grid Gauss-Legcndrc quadrature rule of 
dimension m + n and level p + 1. Following the approach given in Sec. 4, we construct the 
quadrature rule for integration with respect to the joint probability distribution of rj^'^ and C, 
whose nodes and weights we denote as {(^f., Cfe; '^^D' 1 ^ ^ ^ '^^li each iteration as a 
sparse-grid quadrature rule of dimension d + n and level q + 2. This rule is synthesized from 
the family of fully tensorized Gauss-Legendre quadrature rules for integration with respect to 
the probability distribution of ^ and the family of embedded quadratiirc riilcs for integration 
with respect to the probability distribution of tj^'P obtained using the method given in [12]. 

5.9. Selection of the reduced dimension and polynomial degree 

At each iteration, we select the number of terms retained in (67) as the smallest dimension d 
that satisfies 



/ [ \\T(^Ht(:)-T(''''''{tC)f dP^dPc<eiJ [ [ Wr^^PitoW dP^dPc, € N, (71) 

where ei is a prescribed tolerance level. Further, at each iteration, we truncate the chaos 
expansion of the random neutron flux in (70) at the smallest total degree q that satisfies 



\%<^.i{nX)-^^-''-Hv,C)f dPidPi<e2j f [ \\$^'i{r]X)f dP^idPi, V£ € N, (72) 

where €2 is a prescribed tolerance level. Clearly, these criteria may result in the dependence 
of d and q on the number of iterations £. 

5.10. Concluding rem,arks 

Algorithm 1 summarizes an implementation of the problem in which the nonintrusive 
projection method is used for the solution of the subproblems. Although this algorithm uses 
the nonintrusive projection method, we note that the proposed methodology can be readily 
adapted for use with other methods, such as embedded projection and collocation. 

The main feature of the proposed implementation is that it provides a solution of the 
ncutronics subproblcm in a reduced-dimensional space when the reduced chaos expansion with 
random coefficients can extract a low-dimensional representation of the random temperature 
{d < m), while maintaining accuracy. The solution in a reduced-dimensional space can be 
expected to reduce the number of terms required in the chaos expansion of the random 
neutron flux to achieve sufHcient accuracy. Further, it can be expected to reduce the number of 
quadrature nodes required for the nonintrusive projection method to achieve sufficient accuracy 
in the coefficients of the chaos expansion of the random neutron flux. Hence, the solution of the 
neutronics subproblem in a reduced-dimensional space reduces the number of times a sample 
of the neutronics subproblem must be solved, thus in turn lowering the computational cost. 



6. Numerical results 

We obtained numerical results by considering the following parameter values. We assumed 
the reactor to have a length L of 100 [cm]. Further, we assumed a deterministic and position- 
independent heat conductivity k of 100 [J/ (K cms)], ambient temperature T^o of 390 [K], a 
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Input : Error tolerance levels ei and €2] 

Basis {^ipail}0, < |q;| + \(3\ < p} up to total degree p; 
Quadrature rule {(4^, Cfei 'ffe)i 1 < A; < iV} of level p + 1; 



Solve [K + Hi^,)]f' C J = Cfc), iv'-^'^iU, Cfc)) ; 



end 

Compute chaos coordinates of T^^^ using Ti^ = J2k=i^^i^k,<k)vaiik}i'0{CkW; 
end 

dimension reduction 



Compute = Ti^ and Cl^^_ = Ef^iT*""'"^" ^iMTf^^)^; 

Solve the eigenproblem Ef^^lo M^^C^^j,_ W</)^' = ^'jWcf>'/, < |/3| < p - 1; 

Choose d such that ^E^I^ < ^VEf„|+|;3|=o(^i/3)^W-^i0; 



repeat 

measure transformation 

Compute basis {F^V/Si < I7I + |/3| < g} up to total degree q; 
Compute quadrature rule {(tj^, Ct, w^), 1 < A; < i^*} of level q + 2; 
end 

for fc = 1 to do 



Solve [D(T^-'^(^Lci)) +M(T^-^(^l,Cl),Cl)]*^(»7LCD = 
with T^'^'<^(r,t, ct) = T^'^(Ci) + E •=! v^<fc«^''''""'(Ct); 



Compute coordinates of using = Yll=i ^\viXi)^''r{vi}i'0{Ci)wU 

q = q + l; 



repeat 

heat subproblem 

for k = 1 to N do 




end 

neutronics subproblem 

9 = 0; 



end 




end 



until (convergence); 



Algorithm 1: Implementation of the illustration problem. 
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fission energy E[ of 3.0i?-ll [J/neutrons], a fission eross section Sa.rcf of 0.0075 [cm" '^J, a 
neutron-diffusion constant D^ef of 2.2 [cm], a multiplication factor v of 2.2, a neutron source 
strength s of S.Oi^ll [neutrons/ (sem^)]. In addition, the following temperatures were assumed: 
T,ef = 390 [K], T^i„ = 390 [K], and T^^ = 1000 [K]. 

[Figure 2 about here.] 

[Figure 3 about here.] 

Wc used thermal transmittancc and absorption cross section random fields with position- 
independent mean values h = 0.17 [J/(Kcm'^ s)] and Sa.ref = 0.0195 [cm~^], spatial correlation 
lengths ah = 15 [cm] and as = 50 [cm], and coefficients of variation Sh = W% and 5^ = 10%. 
Wc retained m = 10 and n ~ 2 terms in expansions (49) and (50), respectively. 

Figures 2(a) and 3(a) show a few sample paths of the random fields thus obtained. We can 
observe that the samples of the thermal transmittance random field with a/, = 15 [cm] are less 
smooth than those of the absorption cross section random field with as ~ 50 [cm], that is, the 
samples of the former random field exhibit more rapid oscillations with respect to the position 
in the reactor than those of the latter random field. Figures 2(b) and 3(b) show the 10 largest 
eigenvalues of the covariance integral operators whose kernels are given by (51) and (52). We 
can observe that the eigenvalues obtained for the thermal transmittance random field decay 
at a lower rate than those obtained for the absorption cross section random field, indicating 
the adequacy of the truncation of expansions (49) and (50) after m = 10 and n = 2 terms. 

6.1. Monte Carlo sampling implementation 

[Figure 4 about here.] 

First, we carried out a Monte Carlo simulation. We generated MC = 100, 000 pairs of sample 
paths of the thermal transmittance and absorption cross section random fields. Then, for each 
pair of sample paths, we constructed the associated deterministic multiphysics model. We 
solved each model by using the FE method for spatial discretization and Gauss-Seidel iteration. 
We systematically obtained converged results for 40 finite elements and 20 iterations. 

Figure 4 shows a few samples of the random temperature and neutron flux thus obtained. 
The samples of the random temperature (Fig. 4(a)) are smoother than those of the thermal 
transmittance random field (Fig. 2(a)), that is, the former samples exhibit less rapid oscillations 
with respect to the position in the reactor than the latter samples. In [11], we had shown that 
this behavior can be attributed to the large magnitude of the diffusion term of the heat 
subproblem, which reduces the nonuniformity of the samples of the random temperature. 

6.2. PC-based implementation involving dimension reduction and measure transformation 

Next, we implemented the proposed polynomial-chaos-based iterative method involving 
dimension reduction and measure transformation. This implementation corresponded exactly 
to Algorithm 1. We obtained results by setting the total degree p of the chaos expansion of the 
random temperature to 4 and, with reference to (71) and (72), by using a range of values for 
the error tolerance levels ei and 62 to determine the reduced dimension and the total degree of 
the chaos expansion of the random neutron flux at each iteration. We discuss the convergence 
of the results as a function of these error tolerance levels later. Now, we present detailed results 
obtained for ei = 0.01 and €2 = 0.01. 
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[Figure 5 about here.] 

Figure 5 shows the convergence of the iterative method as a function of the number of 

iterations; note that the superscript oo is used in the figure captions to indicate convergence 
with respect to the number of iterations. The iterative method converged hnearly up to 
approximately iteration 1 = 1, after which hnear-solver tolerances became dominant and 
prevented further convergence. All results to follow were obtained at iteration f = 20 and 
can thus be considered to have converged with respect to the number of iterations. 

[Figure 6 about here.] 

[Figure 7 about here.] 

Figure 6 shows a few components of the reduced chaos expansion with random coefficients of 
the random temperature. We can observe that the eigenvalues of the reduced chaos expansion 
with random coefficients of the random temperature (Fig. 6(b)) decay at a higher rate than 
those of the KL decomposition of the thermal transmittance random field (Fig. 2(b)). This 
result is consistent with our earlier observation that the samples of the random temperature 
are smoother than those of the thermal transmittance random field. 

Figure 7 shows the joint and marginal probability density functions of the reduced random 
variables. Clearly, the joint probability density function shows statistical dependence, and the 
marginal probability density functions are not "labeled." 

At iteration = 20, a reduced chaos expansion with random coefficients obtained by 
retaining only d = 2 terms was found to be sufficiently accurate to satisfy (71) for ei = 0.01; 
thus, at this iteration, the measure transformation necessitated the construction of polynomial 
chaos and quadrature rules with respect to the joint probability distribution of the reduced 
random variables rj^'P = {Vi^jvi'^) and the input random variables C = (Ci)C2)- 

[Figure 8 about here.] 

Figure 8 illustrates the proposed computational construction of polynomial chaos with 
respect to ; specifically, the figure shows the polynomial chaos obtained up to a total degree q 
of 2. With reference to Sec. 5.8, the requisite polynomial chaos with respect to x are 
synthesized from the computed polynomial chaos thus obtained with respect to and the 
normalized Legendre polynomials with respect to P(^ by tensorization. 

[Figure 9 about here.] 

Figure 9 illustrates the proposed computational construction of quadrature rules with respect 
to P^; specifically, the figure shows the quadrature rules obtained up to a level A of 4. 
This figure indicates that as A was increased, higher accuracy was required, and thus, a 
quadrature rule with more nodes and weights was systematically obtained. With reference 
to Sec. 5.8, the requisite quadrature rules with respect to P^ x are synthesized from 
the computed quadrature rules obtained with respect to P^ and the fully tensorized Gauss- 
Legendre quadrature rules with respect to P^ by tensorization. 

[Figure 10 about here.] 
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At iteration £ = 20, a chaos expansion truneated at q = 2 was found to be sufficiently 
accurate to satisfy (72) for €2 = 0.01. The representation of the random temperature by a 
chaos expansion of dimension m + n = 12 and total degree p = 4 requires 1, 820 = 16!/12!/4! 
terms; in contrast, the representation of the random neutron flux by a chaos expansion of 
dimension d + n = i and total degree q = 2 requires only 15 = 6!/4!/2! terms. The sparse- 
grid Gauss-Legendre quadrature rule of dimension m + n = 12 and level p + 1 = 5 used to 
compute the coefficients of the chaos expansion of the random temperature has 34, 065 nodes 
and weights; in contrast, the sparse-grid quadrature rule of dimension d + n = 4 and level 
q + 2 = 4: used to compute the coefficients of the chaos expansion of the random neutron flux 
has only 346 nodes and weights. Figure 10 shows a few coefflcients of the solution. 

[Figure 11 about here.] 

Figure 11 shows a few samples of the random temperature and random neutron flux deduced 
from the chaos expansions obtained as the output of the solution algorithm. The samples of 
the input random variables used to synthesize the samples of the random temperature and 
random neutron flux shown in Fig. 11 were identical to those used to generate the samples 
shown in Fig. 4. The similarity of the samples in Figs. 4 and 11 indicates that the surrogate 
model based on polynomial chaos not only provides an accurate global representation of the 
multiphysics model but is also capable of accurately reproducing a sample-wise response. 

6. 3. Sensitivity analysis of the random temperature and random neutron flux 

The chaos expansions in (70) facilitate the following decompositions of the variances of the 
random temperature and random neutron flux: 

= V[ + Vf + V(l^^ , = + V^" + , (73) 

lal + |/3| = l |7|-H/3|=1 
1/3| = 1 1^1 = 1 

^(L)= E WT^^Ww, v^(|,c)= E w^^^Ww- 

l«l + |/3|=l |7| + I/3|=1 

Here, and V^* are the sums of the variances of those terms in the chaos expansions of the 
random temperature and random neutron flux, respectively, which depend only on the input 
random variables ^ that describe the uncertainty in the parameters of the heat subproblem; 
hence, following the approach given in [22- 24], V^-^ and can be interpreted as the portions 
contributed by the uncertainty in the parameters of the heat subproblem to the variances of 
the random temperature and random neutron flux, respectively. Conversely, and V^* are 
the sums of the variances of those terms in the chaos expansions of the random temperature 
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and random neutron flux which depend only on the input random variables ^ that describe 
the uncertainty in the parameters of the neutronics subproblem, and they can be interpreted 
as the portions contributed by the uncertainty in the parameters of the neutronics subproblem 
to the variances of the random temperature and random neutron flux, respectively. Lastly, 
and are the portions contributed by the interaction of the uncertainties in the 

parameters of the heat subproblem and those in the parameters of the neutronics subproblem 
to the variances of the random temperature and random neutron flux, respectively. 

[Table 1 about here.] 

We obtained the values given in Table I; these values indicate that the uncertainty 
in the parameters of the heat subproblem predominantly drives the uncertainty in the 
temperature (55.45%) and that the uncertainty in the parameters of the neutronics subproblem 
predominantly drives the uncertainty in the neutron flux (97.83%). Further, we can observe 
that the uncertainty in the parameters of the heat subproblem is less important in driving 
the uncertainty in the neutron flux than the uncertainty in the parameters of the neutronics 
subproblem is in driving the uncertainty in the temperature (2.15% versus 44.40%). This 
observation can be explained as a consequence of the coupling mechanism adopted in this 
work; with reference to (47) , this coupling mechanism involves the division of the temperature 
by a large reference temperature as it feeds into the coefiicients of the neutronics subproblem, 
thus diminishing the impact that fluctuations in the temperature have on the neutron flux. 

6.4- Convergence analysis 

[Figure 12 about here.] 

[Figure 13 about here.] 

We repeated the simulation based on polynomial chaos for several values of the error 
tolerance levels. Each error tolerance level corresponded to a specific accuracy that the reduced 
chaos expansion with random coefficients of the random temperature and the chaos expansion 
of the random neutron flux were required to maintain at each iteration. Figures 12(a) and 13(a) 
indicate that more terms were retained in these expansions when higher accuracy was required. 

Further, Figs. 12 ((b) and (c)) and 13 ((b) and (c)) indicate that the distance between the 
successive approximations determined by the simulation based on Monte Carlo sampling and 
the simulation based on polynomial chaos remained bounded as the iterations progressed and 
that this distance can be reduced systematically by improving the accuracy of the reduced 
chaos expansion with random coeflacients of the random temperature and the chaos expansion 
of the random neutron flux by decreasing the respective error tolerance levels. 

6.5. Concluding remarks 

The proposed methodology provided the solution of the neutronics subproblem in a reduced- 
dimensional space because the reduced chaos expansion with random coefficients could 
facilitate a low-dimensional representation of the random temperature as it passed from 
the heat subproblem to the neutronics subproblem. While accuracy was maintained, the 
solution in a reduced-dimensional space resulted in computational gains because of the 
following two factors. First, the solution in a reduced-dimensional space facilitated the accurate 
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representation of the random neutron flux by a cliaos expansion that contained only a few 
terms. Second, the coefficients in the chaos expansion of the random neutron flux could be 
computed by using a quadrature rule that had only a few nodes, and thus, the solution of only 
a few samples of the neutronics subproblem was required at each iteration. 



7. Conclusion 

We presented a characterization of information by a reduced chaos expansion with 
random coefficients as this information passes from a subproblem of a stochastic coupled 
problem to another and from iteration to iteration. This expansion provides a reduced- 
dimensional representation of the exchanged information, while maintaining segregation 
between statistically independent sources of uncertainty that stem from different subproblems. 
Further, we presented a measure transformation that allows stochastic expansion methods to 
exploit this dimension reduction to obtain an efficient solution of subproblems in a reduced- 
dimensional space. We showed that owing to the uncertainty source segregation, polynomial 
chaos and quadrature rules required for the implementation of this measure transformation can 
be readily obtained by tensorization. The effectiveness of the methodology was demonstrated 
by considering a multiphysics problem in nuclear engineering. 
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Figure 1. Schematic representation of the problem. 
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Figure 2. Thermal transmittance random field: (a) five samples and (b) ten largest eigenvalues of the 
covariance integral operator whose kernel is given by (51). 
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Figure 3. Absorption cross section random field: (a) five samples and (b) ten largest eigenvalues of the 
covariance integral operator whose kernel is given by (52). 
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Figure 4. Monte Carlo simulation: five samples of the solution. 
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Figure 5. PC-based simulation: convergence with respect to the number of iterations. 
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Figure 6. PC-based simulation: a few components of the reduced chaos expansion with random 
coefficients, T°°'P''^{^, ^) = t°°.p(^) + J2j=i yp^vT'^iO <P''"'-''~\<), of the random temperature. 
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Figure 7. PC-based simulation: probability distribution of the first and second reduced random 
variables of the reduced chaos expansion with random coefficients of the random temperature. 
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Figure 8. PC-based simulation: computed polynomial chaos up to a total degree q of 2. 



Copyright © 2000 John Wiley & Sons, Ltd. 
Prepared using nmeauth.cis 



Int. J. Numer. Meth. Engng 2000; 00:0-0 



FIGURES 



39 



5 
2.5 


-2.5 



-2.5 2.5 

Reduced coordinate r| ^ [-] 



(a) {v'^' , 1 < < 1^^°} for A = 1. 

5r 

2.5- 


-2.5 



-2.5 2.5 

Reduced coordinate r| [-] 



(c) {tj^' , 1 < fc < u-^} for A = 2. 

5r 

2.5- 


-2.5 
-5 



-5 -2.5 2.5 

Reduced coordinate r| ^ [-] 

(e) {T?^'^,l<fc<i/^}for A = 3. 

5| ' ■ ' 

2.5- 





-2.5 



-2.5 2.5 

Reduced coordinate r| [-] 



(g) {r?r' .i<fc<^'nf™^^ 



1 

0.5 


-0.5 











,11 

























10 



20 
Index [ 



30 



40 



(b) 1 < fc < 1^?°} for A = 1. 



1 

0.5 


-0.5 



10 



20 
Index [-] 



30 



40 



(d) ,1 < fc < ly^} for A : 



-1, 











TT.Tf.TT., 























10 20 30 40 

Index [-] 

(f) {w"^'^, 1 < fc < u-^} for A = 3. 

Ir 



0.5 


-0.5 
-1, 























>•••••• 











10 



20 
Index [-] 



30 



40 



(h) ,1 < fc < ly^} for A : 



Figure 9. PC-based simulation: computed quadrature rules up to a level A of 4. 
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Figure 10. PC-based simulation: coefficients of the solution at a; = 10 [cm] . 



Copyright © 2000 John Wiley & Sons, Ltd. 
Prepared using nmeauth.cis 



Int. J. Numer. Meth. Engng 2000; 00:0-0 




20 40 60 80 100 

Position [cm] 

(b) Neutron flux. 
Figure 11. PC-based simulation: five samples of the solution. 
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Figure 12. Convergence analysis: (a) reduced dimension and (b, c) mean-square distance between 
the successive approximations determined by the simulation based on Monte Carlo sampling and 
the simulation based on polynomial chaos for ei = 0.02 (circles), ei = 0.01 (squares), and ei = 
0.005 (diamonds) and £2 = 0.01 as a function of the number of iterations. 
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Figure 13. Convergence analysis: (a) total degree and (b,c) mean-square distance between the 
successive approximations determined by the simulation based on Monte Carlo sampling and the 
simulation based on polynomial chaos for ei = 0.01 and e-z = 0.1 (circles), e-z = 0.01 (squares), 
and 62 = 0.001 (diamonds) as a function of the number of iterations. 
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temperature 


neutron flux 


main etiect 
of uncertainty in data of heat subproblem 


^ = 55.45% 


% = 2.15% 


main effect 

of uncertainty in data of neutronics subproblem 


^ = 44.40% 


% = 97.83% 


interaction 
of uncertainty in data of heat subproblem 
and uncertainty in data of neutronics subproblem 


= 0.02% 


%i - 0.08% 



Table I. Decomposition of the variances of the random temperature and random neutron flux. 
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